home *** CD-ROM | disk | FTP | other *** search
/ Amiga Games: Greatest Hits 1996 / Amiga Games: Greatest Hits 1996.iso / archive / userbox / publicdomain / fractal.lha / fractal / source_original / ENC.C < prev    next >
C/C++ Source or Header  |  1996-07-02  |  41KB  |  946 lines
  1. /**************************************************************************/
  2. /* Encode a byte image using a fractal scheme with a quadtree partition   */
  3. /*                                                                        */
  4. /*       Copyright 1993,1994 Yuval Fisher. All rights reserved.           */
  5. /*                                                                        */
  6. /* Version 0.03 3/14/94                                                   */
  7. /**************************************************************************/
  8.  
  9. #include <stdio.h>
  10. #include <math.h>
  11.  
  12. #define DEBUG 0
  13. #define GREY_LEVELS 255
  14.  
  15. #define bound(a)   ((a) < 0.0 ? 0 : ((a)>255.0? 255 : a))
  16. #define IMAGE_TYPE unsigned char /* may be different in some applications */
  17.  
  18. /* various function declarations to keep compiler warnings away. ANSI     */
  19. /* prototypes can go here, for the hearty.                                */
  20. void fatal();
  21. char *malloc();
  22. char *strcpy(); 
  23.  
  24. /* The following #define allocates an hsize x vsize  matrix of type TYPE */
  25. #define matrix_allocate(matrix, hsize, vsize, TYPE) {\
  26.     TYPE *imptr; \
  27.     int _i; \
  28.     matrix = (TYPE **)malloc((vsize)*sizeof(TYPE *));\
  29.     imptr = (TYPE*)malloc((long)(hsize)*(long)(vsize)*sizeof(TYPE));\
  30.     if (imptr == NULL) \
  31.         fatal("\nNo memory in matrix allocate."); \
  32.     for (_i = 0; _i<vsize; ++_i, imptr += hsize) \
  33.         matrix[_i] = imptr; \
  34.  }
  35.  
  36. #define swap(a,b,TYPE)           {TYPE _temp; _temp=b; b=a; a= _temp;}
  37.  
  38. IMAGE_TYPE **image;          /* The input image data                      */ 
  39. double     **domimage[4];    /* Decimated input image used for domains    */
  40.  
  41. double max_scale = 1.0;      /* Maximum allowable grey level scale factor */
  42.  
  43. int    s_bits = 5,           /* Number of bits used to store scale factor */
  44.        o_bits = 7,           /* Number of bits used to store offset       */
  45.        min_part = 4,         /* Min and max _part determine a range of    */
  46.        max_part = 6,         /* Range sizes from hsize>>min to hsize>>max */
  47.        dom_step = 1,         /* Density of domains relative to size       */
  48.        dom_step_type = 0,    /* Flag for dom_step a multiplier or divisor */
  49.        dom_type = 0,         /* Method of generating domain pool 0,1,2..  */
  50.        only_positive = 0,    /* A flag specifying use of positive scaling */
  51.        subclass_search = 0,  /* A flag specifying classes searched        */
  52.        fullclass_search = 0, /* A flag specifying classes searched        */
  53.        *bits_needed,         /* Number of bits to encode domain position. */
  54.        zero_ialpha,          /* The const ialpha when alpha = 0           */
  55.        max_exponent;         /* The max power of 2 side of square image   */
  56.                              /* that fits in our input image.             */
  57.  
  58.  
  59.                              /* The class_transform gives the transforms  */
  60.                              /* between classification numbers for        */
  61.                              /* negative scaling values, when brightest   */
  62.                              /* becomes darkest, etc...                   */
  63. int    class_transform[2][24] = {23,17,21,11,15,9,22,16,19,5,13,3,20,10,18,
  64.                                  4,7,1,14,8,12,2,6,0,
  65.                                  16,22,10,20,8,14,17,23,4,18,2,12,11,21,5,
  66.                                  19,0,6,9,15,3,13,1,7};
  67.  
  68.                              /* rot_transform gives the rotations for     */
  69.                              /* domains with negative scalings.           */
  70. int    rot_transform[2][8] = {7,4,5,6,1,2,3,0, 2,3,0,1,6,7,4,5};
  71.  
  72. struct domain_data {
  73.        int *no_h_domains,    /* The number of domains horizontally for    */
  74.            *no_v_domains,    /* each size.                                */
  75.            *domain_hsize,    /* The size of the domain.                   */
  76.            *domain_vsize,    /* The size of the domain.                   */
  77.            *domain_hstep,    /* The density of the domains.               */
  78.            *domain_vstep;    /* The density of the domains.               */
  79. struct domain_pixels {       /* This is a three (sigh) index array that   */
  80.        int dom_x, dom_y;     /* dynamically allocated. The first index is */
  81.        double sum,sum2;      /* the domain size, the other are two its    */
  82.        int sym;              /* position. It contains the sum and sum^2   */
  83. } ***pixel;                  /* of the pixel values in the domains, which */
  84. } domain;                    /* are computed just once.                   */ 
  85.                              
  86.  
  87. struct classified_domain {             /* This is a list which containes  */
  88.        struct domain_pixels *the;      /* pointers to  the domain data    */
  89.        struct classified_domain *next; /* in the structure above. There   */
  90. } **the_domain[3][24];                 /* are three classes with 24 sub-  */
  91.                                        /* classes. Using this array, only */
  92.                                        /* domains and ranges in the same  */
  93.                                        /* class are compared..            */
  94.                                        /* The first pointer points to the */
  95.                                        /* domain size the the second to   */
  96.                                        /* list of domains.                */
  97.  
  98. FILE *output;                /* Output FILE containing compressed data    */
  99.  
  100. main(argc,argv)
  101. int argc;
  102. char **argv;
  103. /* Usage: quadfrac [tol [inputfilename [outputfilename [hsize [vsize]]]]]  */
  104. {
  105.     /* Defaults are set initially */
  106.     double         tol = 8.0;            /* Tolerance value for quadtree.  */
  107.     char           inputfilename[200];
  108.     char           outputfilename[200];
  109.     int            i,j,k,
  110.                    hsize = -1,           /* The horizontal and vertical    */
  111.                    vsize = -1;           /* size of the input image.       */
  112.     long           stripchar=0;          /* chars to ignore in input file. */
  113.     FILE           *input; 
  114.  
  115.     inputfilename[0] = 1;  /* We initially set the input to this and */
  116.     outputfilename[0] = 1; /* then check if the input/output names   */
  117.                            /* have been set below.                   */
  118.  
  119.     /* scan through the input line and read in the arguments */
  120.     for (i=1; i<argc; ++i) 
  121.        if (argv[i][0] != '-' ) 
  122.             if (inputfilename[0] == 1)
  123.                 strcpy(inputfilename, argv[i]);
  124.             else if (outputfilename[0] == 1)
  125.                 strcpy(outputfilename, argv[i]);
  126.             else;
  127.        else { /* we have a flag */
  128.             if (strlen(argv[i]) == 1) break;
  129.             switch(argv[i][1]) {
  130.                case 't': tol = atof(argv[++i]);
  131.                          break;
  132.                case 'S': stripchar = atoi(argv[++i]);
  133.                          break;
  134.                case 'x': 
  135.                case 'w': hsize = atoi(argv[++i]);
  136.                          break;
  137.                case 'y': 
  138.                case 'h': vsize = atoi(argv[++i]);
  139.                          break;
  140.                case 'D': dom_type = atoi(argv[++i]);
  141.                          break;
  142.                case 'd': dom_step = atoi(argv[++i]);
  143.                          if (dom_step < 0 || dom_step > 15) 
  144.                             fatal("\n Bad domain step.");
  145.                          break;
  146.                case 's': s_bits = atoi(argv[++i]);
  147.                          break;
  148.                case 'o': o_bits = atoi(argv[++i]);
  149.                          break;
  150.                case 'm': min_part = atoi(argv[++i]);
  151.                          break;
  152.                case 'M': max_part = atoi(argv[++i]);
  153.                          break;
  154.                case 'e': dom_step_type= 1;
  155.                          break;
  156.                case 'p': only_positive = 1;
  157.                          break;
  158.                case 'f': subclass_search = 1;
  159.                          break;
  160.                case 'F': fullclass_search = 1;
  161.                          break;
  162.                case 'N': max_scale = atof(argv[++i]);
  163.                          break;
  164.                case '?':
  165.                case 'H':
  166.                default:
  167.            printf("\nUsage: enc -[options] [inputfile [outputfile]]");
  168.            printf("\nOptions are: (# = number), defaults show in ()");
  169.            printf("\n  -t # tolerance criterion for fidelity. (%lf)", tol);
  170.            printf("\n  -m # minimum quadtree partitions.  (%d)",min_part);
  171.            printf("\n  -M # maximum quadtree partitions. (%d)",max_part);
  172.            printf("\n  -S # number of input bytes to ignore. (%ld)",stripchar);
  173.            printf("\n  -w # width (horizontal size) of input data. (256)");
  174.            printf("\n  -h # height (vertical size) of input data. (256)");
  175.            printf("\n  -d # domain step size. (%d)", dom_step);
  176.            printf("\n  -D # method {0,1,2} for domain pool (%d)",dom_type);
  177.            printf("\n  -s # number of scaling quantizing bits. (%d)",s_bits);
  178.            printf("\n  -o # number of offset quantizing bits. (%d)",o_bits);
  179.            printf("\n  -N # maximum scaling in encoding. (%lf)",max_scale);
  180.            printf("\n  -e   domain step as multiplier not divisor. (off)");
  181.            printf("\n  -p   use only positive scaling (for speed). (off)");
  182.            printf("\n  -f   search 24 domain classes (for fidelity). (off)");
  183.            printf("\n  -F   search 3 domain classes (for fidelity). (off)");
  184.            fatal("\n       -F and -f can be used together.");
  185.             }
  186.        }
  187.    
  188.     if (inputfilename[0] == 1) strcpy(inputfilename, "lenna.dat");
  189.     if (outputfilename[0] == 1) strcpy(outputfilename, "lenna.trn");
  190.  
  191.     if (hsize == -1) 
  192.         if (vsize == -1) hsize = vsize = 256;
  193.         else hsize = vsize;
  194.     else 
  195.         if (vsize == -1) vsize = hsize;
  196.  
  197.     /* allocate memory for the input image. Allocating one chunck saves  */
  198.     /* work and time later.                                              */
  199.     matrix_allocate(image, hsize, vsize, IMAGE_TYPE)
  200.     matrix_allocate(domimage[0], hsize/2, vsize/2, double)
  201.     matrix_allocate(domimage[1], hsize/2, vsize/2, double)
  202.     matrix_allocate(domimage[2], hsize/2, vsize/2, double)
  203.     matrix_allocate(domimage[3], hsize/2, vsize/2, double)
  204.  
  205.     /* max_ & min_ part are variable, so this must be run time allocated */
  206.     bits_needed = (int *)malloc(sizeof(int)*(1+max_part-min_part));
  207.  
  208.     if ((input = fopen(inputfilename, "r")) == NULL)
  209.         fatal("Can't open input file.");
  210.  
  211.     /* skip the first  stripchar  chars */ 
  212.     fseek(input, stripchar, 0); 
  213.     i = fread(image[0], sizeof(IMAGE_TYPE), hsize*vsize, input);
  214.     fclose(input);
  215.  
  216.     if (i < hsize*vsize) 
  217.         fatal("Not enough input data in the input file.");
  218.     else
  219.         printf("%dx%d=%d pixels read from %s.", hsize,vsize,i,inputfilename);
  220.  
  221.     /* allcate memory for domain data and initialize it */
  222.     compute_sums(hsize,vsize);
  223.  
  224.     if ((output = fopen(outputfilename, "w")) == NULL)
  225.         fatal("Can't open output file.");
  226.  
  227.     /* output some data into the outputfile.                       */
  228.     pack(4,(long)min_part,output);
  229.     pack(4,(long)max_part,output);
  230.     pack(4,(long)dom_step,output);
  231.     pack(1,(long)dom_step_type,output);
  232.     pack(2,(long)dom_type,output);
  233.     pack(12,(long)hsize,output);
  234.     pack(12,(long)vsize,output);
  235.  
  236.     /* This is the quantized value of zero scaling.. needed later */
  237.     zero_ialpha = 0.5 + (max_scale)/(2.0*max_scale)*(1<<s_bits);
  238.  
  239.     /* The following routine takes a rectangular image and calls the */
  240.     /* quadtree routine to encode square sum-images in it.           */
  241.     /* the tolerance is a parameter since in some applications different */
  242.     /* regions of the image may need to be compressed to different tol's */
  243.     printf("\nEncoding Image.....");
  244.     fflush(stdout);
  245.     partition_image(0, 0, hsize,vsize, tol);
  246.     printf("Done.");
  247.     fflush(stdout);
  248.  
  249.     /* stuff the last byte if needed */
  250.     pack(-1,(long)0,output);
  251.     
  252.     fclose(output);
  253.     i = pack(-2,(long)0,output);
  254.     printf("\n Compression = %lf from %d bytes written in %s.\n",
  255.             (double)(hsize*vsize)/(double)i, i, outputfilename);
  256.  
  257.     /* Free allocated memory*/
  258.     free(bits_needed);
  259.     free(domimage[0]);
  260.     free(domimage[1]);
  261.     free(domimage[2]);
  262.     free(domimage[3]);
  263.     free(domain.no_h_domains);
  264.     free(domain.no_v_domains);
  265.     free(domain.domain_hsize);
  266.     free(domain.domain_vsize);
  267.     free(domain.domain_hstep);
  268.     free(domain.domain_vstep);
  269.     for (i=0; i <= max_part-min_part; ++i) 
  270.         free(domain.pixel[i]);
  271.     free(domain.pixel);
  272.     free(image[0]);
  273.     for (i=0; i <= max_part-min_part; ++i) 
  274.     for (k=0; k<3; ++k) 
  275.      for (j=0; j<24; ++j) list_free(the_domain[k][j][i]);
  276.     return(0);
  277. }
  278.  
  279. /* ************************************************************** */
  280. /* free memory allocated in the list structure the_domain         */
  281. /* ************************************************************** */
  282. list_free(node)
  283. struct classified_domain *node;
  284. {
  285.      if (node->next != NULL)
  286.        list_free(node->next);
  287.      free(node);
  288. }
  289.  
  290. /* ************************************************************** */
  291. /* return the average pixel value of a region of the image.       */
  292. /* ************************************************************** */
  293. void average(x,y,xsize,ysize, psum, psum2)
  294. int x,y,xsize,ysize;
  295. double *psum, *psum2;
  296. {
  297.     register int i,j,k;
  298.     register double pixel;
  299.     *psum = *psum2 = 0.0;
  300.     k = ((x%2)<<1) + y%2;
  301.     x >>= 1; y >>= 1;
  302.     xsize >>= 1; ysize >>= 1;
  303.     for (i=x; i<x+xsize; ++i)
  304.     for (j=y; j<y+ysize; ++j) {
  305.        pixel = domimage[k][j][i];
  306.        *psum += pixel;
  307.        *psum2 += pixel*pixel;
  308.     }
  309. }
  310.  
  311. /* ************************************************************** */
  312. /* return the average pixel value of a region of the image. This  */
  313. /* routine differs from the previous in one slight way. It does   */
  314. /* not average 2x2 sub-images to pixels. This is needed for clas- */
  315. /* sifying ranges rather than domain where decimation is needed.  */
  316. /* ************************************************************** */
  317. void average1(x,y,xsize,ysize, psum, psum2)
  318. int x,y,xsize,ysize;
  319. double *psum, *psum2;
  320. {
  321.     register int i,j;
  322.     register double pixel;
  323.     *psum = *psum2 = 0.0;
  324.  
  325.     for (i=x; i<x+xsize; ++i)
  326.     for (j=y; j<y+ysize; ++j) {
  327.        pixel = (double)image[j][i];
  328.        *psum += pixel;
  329.        *psum2 += pixel*pixel;
  330.     }
  331. }
  332.  
  333. /* ************************************************************** */
  334. /* Take a region of the image at x,y and classify it.             */
  335. /* The four quadrants of the region are ordered from brightest to */
  336. /* least bright average value, then it is rotated into one of the */
  337. /* three cannonical orientations possible with the brightest quad */
  338. /* in the upper left corner.                                      */
  339. /* The routine returns two indices that are class numbers: pfirst */
  340. /* and psecond; the symmetry operation that bring the square into */
  341. /* cannonical position; and the sum and sum^2 of the pixel values */
  342. /* ************************************************************** */
  343. classify(x, y, xsize, ysize, pfirst, psecond, psym, psum, psum2, type)
  344. int x,y,xsize,ysize,   /* position, size of subimage to be classified */
  345.     *pfirst, *psecond, /* returned first and second class numbers     */
  346.     *psym;             /* returned symmetry operation that brings the */
  347.                        /* subimage to cannonical position.            */
  348. double *psum, *psum2;  /* returned sum and sum^2 of pixel values      */
  349. int type;              /* flag for decimating (for domains) or not    */ 
  350. {
  351.  
  352.     int order[4], i,j; 
  353.     double a[4],a2[4];
  354.     void (*average_func)();
  355.     
  356.     if (type == 2) average_func = average; else average_func = average1;
  357.  
  358.     /* get the average values of each quadrant                         */
  359.  
  360.  
  361.     (*average_func)(x,y,                 xsize/2,ysize/2,   &a[0], &a2[0]);
  362.     (*average_func)(x,y+ysize/2,         xsize/2,ysize/2,   &a[1], &a2[1]);
  363.     (*average_func)(x+xsize/2,y+ysize/2, xsize/2,ysize/2,   &a[2], &a2[2]);
  364.     (*average_func)(x+xsize/2,y,         xsize/2,ysize/2,   &a[3], &a2[3]);
  365.  
  366.     *psum = a[0] + a[1] + a[2] + a[3];
  367.     *psum2 =  a2[0] + a2[1] + a2[2] + a2[3];
  368.  
  369.     for (i=0; i<4; ++i) {
  370.          /* after the sorting below order[i] is the i-th brightest   */
  371.          /* quadrant.                                                */
  372.          order[i] = i; 
  373.          /* convert a2[] to store the variance of each quadrant      */
  374.          a2[i] -= (double)(1<<(2*type))*a[i]*a[i]/(double)(xsize*ysize);
  375.     }
  376.  
  377.     /* Now order the average value and also in order[],   which will */
  378.     /* then tell us the indices (in a[]) of the brightest to darkest */
  379.     for (i=2; i>=0; --i)
  380.     for (j=0; j<=i; ++j)
  381.         if (a[j]<a[j+1]) {
  382.             swap(order[j], order[j+1],int)
  383.             swap(a[j], a[j+1],double)
  384.     }
  385.  
  386.     /* because of the way we ordered the a[] the rotation can be */
  387.     /* read right off of order[]. That will make the brightest   */
  388.     /* quadrant be in the upper left corner. But we must still   */ 
  389.     /* decide which cannonical class the image portion belogs    */
  390.     /* to and whether to do a flip or just a rotation. This is   */
  391.     /* the following table summarizes the horrid lines below     */
  392.     /* order      class            do a rotation                 */
  393.     /* 0,2,1,3      0                   0                        */
  394.     /* 0,2,3,1      0                   1                        */
  395.     /* 0,1,2,3      1                   0                        */
  396.     /* 0,3,2,1      1                   1                        */
  397.     /* 0,1,3,2      2                   0                        */
  398.     /* 0,3,1,2      2                   1                        */
  399.  
  400.     *psym = order[0]; 
  401.     /* rotate the values */
  402.     for (i=0; i<4; ++i)
  403.         order[i] = (order[i] - (*psym) + 4)%4; 
  404.  
  405.     for (i=0; order[i] != 2; ++i); 
  406.     *pfirst = i-1;
  407.     if (order[3] == 1 || (*pfirst == 2 && order[2] == 1)) *psym += 4;
  408.  
  409.     /* Now further classify the sub-image by the variance of its    */
  410.     /* quadrants. This give 24 subclasses for each of the 3 classes */
  411.     for (i=0; i<4; ++i) order[i] = i; 
  412.  
  413.     for (i=2; i>=0; --i)
  414.     for (j=0; j<=i; ++j)
  415.         if (a2[j]<a2[j+1]) {
  416.             swap(order[j], order[j+1],int)
  417.             swap(a2[j], a2[j+1],double)
  418.     }
  419.  
  420.     /* Now do the symmetry operation */
  421.     for (i=0; i<4; ++i)
  422.         order[i] = (order[i] - (*psym%4) + 4)%4; 
  423.     if (*psym > 3) 
  424.         for (i=0; i<4; ++i)
  425.            if (order[i]%2) order[i] = (2 + order[i])%4;
  426.  
  427.     /* We want to return a class number from 0 to 23 depending on */
  428.     /* the ordering of the quadrants according to their variance  */
  429.     *psecond = 0;
  430.     for (i=2; i>=0; --i)
  431.     for (j=0; j<=i; ++j)
  432.         if (order[j] > order[j+1]) {
  433.             swap(order[j],order[j+1], int);
  434.             if (order[j] == 0 || order [j+1] == 0) 
  435.                  *psecond += 6;
  436.             else if (order[j] == 1 || order [j+1] == 1) 
  437.                  *psecond += 2;
  438.             else if (order[j] == 2 || order [j+1] == 2) 
  439.                  *psecond += 1;
  440.         }
  441. }
  442.  
  443. /* ************************************************************ */
  444. /* Compute sum and sum^2 of pixel values in domains for use in  */
  445. /* the rms computation later. Since a domain is compared with   */
  446. /* many ranges, doing this just once saves a lot of computation */
  447. /* This routine also fills a list structure with the domains    */
  448. /* as they are classified and creates the memory for the domain */
  449. /* data in a matrix.                                            */
  450. /* ************************************************************ */
  451. compute_sums(hsize,vsize)
  452. int hsize,vsize;
  453. {
  454.     int i,j,k,l,
  455.          domain_x, 
  456.          domain_y, 
  457.          first_class, 
  458.          second_class,
  459.          domain_size, 
  460.          domain_step_size, 
  461.          size,  
  462.          x_exponent, 
  463.          y_exponent;
  464.  
  465.     struct classified_domain *node;
  466.  
  467.     printf("\nComputing domain sums... ");
  468.     fflush(stdout); 
  469.  
  470.     /* pre-decimate the image into domimage to avoid having to      */
  471.     /* do repeated averaging of 2x2 pixel groups.                   */
  472.     /* There are 4 ways to decimate the image, depending on the     */
  473.     /* location of the domain, odd or even address.                 */
  474.     for (i=0; i<2; ++i)
  475.     for (j=0; j<2; ++j)
  476.     for (k=i; k<hsize-i; k += 2)
  477.     for (l=j; l<vsize-j; l += 2) 
  478.         domimage[(i<<1)+j][l>>1][k>>1] = 
  479.                      ((double)image[l][k] + (double)image[l+1][k+1] +
  480.                      (double)image[l][k+1] + (double)image[l+1][k])*0.25;
  481.  
  482.  
  483.     /* Allocate memory for the sum and sum^2 of domain pixels        */
  484.     /* We first compute the size of the largest square that fits in  */
  485.     /* the image.                                                    */
  486.     x_exponent = (int)floor(log((double)hsize)/log(2.0));
  487.     y_exponent = (int)floor(log((double)vsize)/log(2.0));
  488.    
  489.     /* exponent is min of x_ and y_ exponent */
  490.     max_exponent = (x_exponent > y_exponent ? y_exponent : x_exponent);
  491.  
  492.     /* size is the size of the largest square that fits in the image */
  493.     /* It is used to compute the domain and range sizes.             */
  494.     size =  1<<max_exponent; 
  495.  
  496.     if (max_exponent < max_part)
  497.       fatal("Reduce maximum number of quadtree partitions.");
  498.     if (max_exponent-2 < max_part)
  499.       printf("\nWarning: so many quadtree partitions yield absurd ranges.");
  500.  
  501.     i = max_part - min_part + 1;
  502.     domain.no_h_domains = (int *)malloc(sizeof(int)*i);
  503.     domain.no_v_domains = (int *)malloc(sizeof(int)*i);
  504.     domain.domain_hsize = (int *)malloc(sizeof(int)*i);
  505.     domain.domain_vsize = (int *)malloc(sizeof(int)*i);
  506.     domain.domain_hstep = (int *)malloc(sizeof(int)*i);
  507.     domain.domain_vstep = (int *)malloc(sizeof(int)*i);
  508.  
  509.     domain.pixel= (struct domain_pixels ***)
  510.               malloc(i*sizeof(struct domain_pixels **));
  511.     if (domain.pixel == NULL) fatal("No memory for domain pixel sums.");
  512.  
  513.     for (i=0; i <= max_part-min_part; ++i) {
  514.        /* first compute how many domains there are horizontally */
  515.        domain.domain_hsize[i] = size >> (min_part+i-1);
  516.        if (dom_type == 2) 
  517.              domain.domain_hstep[i] = dom_step; 
  518.        else if (dom_type == 1)
  519.              if (dom_step_type == 1) 
  520.                domain.domain_hstep[i] = (size >> (max_part - i-1))*dom_step;
  521.              else 
  522.                domain.domain_hstep[i] = (size >> (max_part - i-1))/dom_step;
  523.        else 
  524.              if (dom_step_type == 1) 
  525.                domain.domain_hstep[i] = domain.domain_hsize[i]*dom_step;
  526.              else 
  527.                domain.domain_hstep[i] = domain.domain_hsize[i]/dom_step;
  528.  
  529.        domain.no_h_domains[i] = 1+(hsize-domain.domain_hsize[i])/
  530.                                                   domain.domain_hstep[i];
  531.        /* bits_needed[i][0] = ceil(log(domain.no_h_domains[i])/log(2.0));  */
  532.  
  533.        /* now compute how many domains there are vertically. The sizes  */
  534.        /* are the same for square domains, but not for rectangular ones */
  535.        domain.domain_vsize[i] = size >> (min_part+i-1);
  536.        if (dom_type == 2) 
  537.              domain.domain_vstep[i] = dom_step; 
  538.        else if (dom_type == 1)
  539.              if (dom_step_type == 1) 
  540.                domain.domain_vstep[i] = (size >> (max_part - i-1))*dom_step;
  541.              else
  542.                domain.domain_vstep[i] = (size >> (max_part - i-1))/dom_step;
  543.        else 
  544.              if (dom_step_type == 1) 
  545.                domain.domain_vstep[i] = domain.domain_vsize[i]*dom_step;
  546.              else
  547.                domain.domain_vstep[i] = domain.domain_vsize[i]/dom_step;
  548.  
  549.        domain.no_v_domains[i] = 1+(vsize-domain.domain_vsize[i])/
  550.                                                   domain.domain_vstep[i];
  551.  
  552.        /* Now compute the number of bits needed to store the domain data */
  553.        bits_needed[i] = ceil(log((double)domain.no_h_domains[i]*
  554.                                  (double)domain.no_v_domains[i])/log(2.0)); 
  555.  
  556.        matrix_allocate(domain.pixel[i], domain.no_h_domains[i],
  557.                      domain.no_v_domains[i], struct domain_pixels) 
  558.  
  559.     }
  560.  
  561.     /* allocate and zero the list containing the classified domain data */
  562.     i = max_part - min_part + 1;
  563.     for (first_class = 0; first_class < 3; ++first_class)
  564.     for (second_class = 0; second_class < 24; ++second_class) {
  565.        the_domain[first_class][second_class] = 
  566.                            (struct classified_domain **) 
  567.                            malloc(i*sizeof(struct classified_domain *));
  568.        for (j=0; j<i; ++j)
  569.                     the_domain[first_class][second_class][j] = NULL;
  570.     }
  571.  
  572.     /* Precompute sum and sum of squares for domains                 */
  573.     /* This part can get faster for overlapping domains if repeated  */
  574.     /* sums are avoided                                              */
  575.     for (i=0; i <= max_part-min_part; ++i) {
  576.       for (j=0,domain_x=0; j<domain.no_h_domains[i]; ++j,
  577.              domain_x+=domain.domain_hstep[i]) 
  578.       for (k=0,domain_y=0; k<domain.no_v_domains[i]; ++k,
  579.              domain_y+=domain.domain_vstep[i]) {
  580.         classify(domain_x, domain_y, 
  581.                  domain.domain_hsize[i], 
  582.                  domain.domain_vsize[i], 
  583.                      &first_class, &second_class,
  584.                      &domain.pixel[i][k][j].sym, 
  585.                      &domain.pixel[i][k][j].sum, 
  586.                      &domain.pixel[i][k][j].sum2, 2);
  587.  
  588.         /* When the domain data is referenced from the list, we need to */
  589.         /* know where the domain is.. so we have to store the position  */
  590.         domain.pixel[i][k][j].dom_x = j;
  591.         domain.pixel[i][k][j].dom_y = k;
  592.         node = (struct classified_domain *)
  593.                                 malloc(sizeof(struct classified_domain));
  594.         
  595.         /* put this domain in the classified list structure */
  596.         node->the = &domain.pixel[i][k][j];
  597.         node->next = the_domain[first_class][second_class][i];
  598.         the_domain[first_class][second_class][i] = node;
  599.       }
  600.     }
  601.  
  602.     /* Now we make sure no domain class is actually empty.  */
  603.     for (i=0; i <= max_part-min_part; ++i) 
  604.     for (first_class = 0; first_class < 3; ++first_class)
  605.     for (second_class = 0; second_class < 24; ++second_class) 
  606.        if (the_domain[first_class][second_class][i] == NULL) {
  607.            node = (struct classified_domain *)
  608.                                 malloc(sizeof(struct classified_domain));
  609.            node->the = &domain.pixel[i][0][0];
  610.            node->next = NULL;
  611.            the_domain[first_class][second_class][i] = node;
  612.        }
  613.  
  614.     printf("Done.");
  615.     fflush(stdout); 
  616. }
  617.  
  618. /* ************************************************************ */
  619. /* pack value using size bits and output into foutf */
  620. /* ************************************************************ */
  621. int pack(size, value, foutf)
  622. int size; long int value;
  623. FILE *foutf;
  624. {
  625.      int i;
  626.      static int ptr = 1, /* how many bits are packed in sum so far */
  627.                 sum = 0, /* packed bits */ 
  628.                 num_of_packed_bytes = 0; /* total bytes written out */
  629.  
  630.     /* size == -1 means we are at the end, so write out what is left */
  631.     if (size == -1 && ptr != 1) {
  632.         fputc(sum<<(8-ptr), foutf);
  633.         ++num_of_packed_bytes;
  634.         return(0);
  635.     }
  636.     
  637.     /* size == -2 means we want to know how many bytes we have written */
  638.     if (size == -2) 
  639.             return(num_of_packed_bytes);
  640.     
  641.     for (i=0; i<size; ++i, ++ptr, value = value>>1, sum = sum<<1) {
  642.         if (value & 1) sum |= 1;
  643.  
  644.          if (ptr == 8) {
  645.             fputc(sum,foutf);
  646.             ++num_of_packed_bytes;
  647.             sum=0;
  648.             ptr=0;
  649.         }
  650.      }
  651. }
  652.  
  653. /* ************************************************************ */
  654. /* Compare a range to a domain and return rms and the quantized */
  655. /* scaling and offset values (pialhpa, pibeta).                 */
  656. /* ************************************************************ */
  657. double compare(atx,aty, xsize, ysize, depth, rsum, rsum2, dom_x,dom_y, 
  658.                                   sym_op, pialpha,pibeta)
  659. int atx, aty, xsize, ysize, depth, dom_x, dom_y, sym_op, *pialpha, *pibeta;
  660. double rsum, rsum2;
  661. {
  662.     int i, j, i1, j1, k,
  663.         domain_x, 
  664.         domain_y;        /* The domain position                */
  665.  
  666.     double pixel,
  667.            det,          /* determinant of solution */
  668.            dsum,         /* sum of domain values */
  669.            rdsum = 0,    /* sum of range*domain values   */
  670.            dsum2,        /* sum of domain^2 values   */
  671.            w2 = 0,       /* total number of values tested */
  672.            rms = 0,      /* root means square difference */
  673.            alpha,        /* the scale factor */
  674.            beta;         /* The offset */
  675.  
  676.  
  677.     
  678.     /* xsize = hsize >> depth; */
  679.     /* ysize = vsize >> depth; */
  680.     w2 = xsize * ysize;
  681.  
  682.     dsum = domain.pixel[depth-min_part][dom_y][dom_x].sum;
  683.     dsum2 = domain.pixel[depth-min_part][dom_y][dom_x].sum2;
  684.     domain_x = (dom_x * domain.domain_hstep[depth-min_part]);
  685.     domain_y = (dom_y * domain.domain_vstep[depth-min_part]);
  686.     k = ((domain_x%2)<<1) + domain_y%2;
  687.     domain_x >>= 1;
  688.     domain_y >>= 1;
  689.  
  690.     /* The main statement in this routine is a switch statement which */
  691.     /* scans through the domain and range to compare them. The loop   */
  692.     /* center is the same so we #define it for easy modification      */
  693. #define COMPUTE_LOOP              {                                 \
  694.         pixel = domimage[k][j1][i1];                                \
  695.         rdsum += image[j][i]*pixel;                                 \
  696.         }
  697.     
  698.      switch(sym_op) {
  699.          case 0: for (i=atx, i1 = domain_x; i<atx+xsize; ++i, ++i1)
  700.                  for (j=aty, j1 = domain_y; j<aty+ysize; ++j, ++j1) 
  701.                      COMPUTE_LOOP
  702.                  break;
  703.          case 1: for (j=aty, i1 = domain_x; j<aty+ysize; ++j, ++i1)
  704.                  for (i=atx+xsize-1, j1 = domain_y; i>=atx; --i, ++j1) 
  705.                      COMPUTE_LOOP
  706.                  break;
  707.          case 2: for (i=atx+xsize-1, i1 = domain_x; i>=atx; --i, ++i1)
  708.                  for (j=aty+ysize-1, j1 = domain_y; j>=aty; --j, ++j1) 
  709.                      COMPUTE_LOOP
  710.                  break;
  711.          case 3: for (j=aty+ysize-1, i1 = domain_x; j>=aty; --j, ++i1)
  712.                  for (i=atx, j1 = domain_y; i<atx+xsize; ++i, ++j1) 
  713.                      COMPUTE_LOOP
  714.                  break;
  715.          case 4: for (j=aty, i1 = domain_x; j<aty+ysize; ++j, ++i1)
  716.                  for (i=atx, j1 = domain_y; i<atx+xsize; ++i, ++j1) 
  717.                      COMPUTE_LOOP
  718.                  break;
  719.          case 5: for (i=atx, i1 = domain_x; i<atx+xsize; ++i, ++i1)
  720.                  for (j=aty+ysize-1, j1 = domain_y; j>=aty; --j, ++j1) 
  721.                      COMPUTE_LOOP
  722.                  break;
  723.          case 6: for (j=aty+ysize-1, i1 = domain_x; j>=aty; --j, ++i1)
  724.                  for (i=atx+xsize-1, j1 = domain_y; i>=atx; --i, ++j1) 
  725.                      COMPUTE_LOOP
  726.                  break;
  727.          case 7: for (i=atx+xsize-1, i1 = domain_x; i>=atx; --i, ++i1)
  728.                  for (j=aty, j1 = domain_y; j<aty+ysize; ++j, ++j1) 
  729.                      COMPUTE_LOOP
  730.                  break;
  731.      }
  732.     
  733.      det = (xsize*ysize)*dsum2 - dsum*dsum;
  734.     
  735.      if (det == 0.0) 
  736.          alpha = 0.0; 
  737.       else 
  738.          alpha = (w2*rdsum - rsum*dsum)/det;
  739.         
  740.      if  (only_positive && alpha < 0.0) alpha = 0.0;
  741.      *pialpha = 0.5 + (alpha + max_scale)/(2.0*max_scale)*(1<<s_bits);
  742.      if (*pialpha < 0) *pialpha = 0;
  743.      if (*pialpha >= 1<<s_bits) *pialpha = (1<<s_bits)-1;
  744.  
  745.      /* Now recompute alpha back */
  746.      alpha = (double)*pialpha/(double)(1<<s_bits)*(2.0*max_scale)-max_scale;
  747.         
  748.      /* compute the offset */
  749.      beta = (rsum - alpha*dsum)/w2;
  750.  
  751.      /* Convert beta to an integer */
  752.      /* we use the sign information of alpha to pack efficiently */
  753.      if (alpha > 0.0)  beta += alpha*GREY_LEVELS;
  754.      *pibeta = 0.5 + beta/
  755.              ((1.0+fabs(alpha))*GREY_LEVELS)*((1<<o_bits)-1);
  756.      if (*pibeta< 0) *pibeta = 0;
  757.      if (*pibeta>= 1<<o_bits) *pibeta = (1<<o_bits)-1;
  758.  
  759.      /* Recompute beta from the integer */
  760.      beta = (double)*pibeta/(double)((1<<o_bits)-1)*
  761.                ((1.0+fabs(alpha))*GREY_LEVELS);
  762.      if (alpha > 0.0) beta  -= alpha*GREY_LEVELS;
  763.  
  764.      /* Compute the rms based on the quantized alpha and beta! */
  765.      rms= sqrt((rsum2 + alpha*(alpha*dsum2 - 2.0*rdsum + 2.0*beta*dsum) +
  766.                  beta*(beta*w2 - 2.0*rsum))/w2);
  767.     
  768.      return(rms);
  769. }
  770.  
  771. /* ************************************************************ */
  772. /* Recursively partition an image, computing the best transfoms */
  773. /* ************************************************************ */
  774. quadtree(atx,aty,xsize,ysize,tol,depth)
  775. int atx, aty, xsize, ysize, depth; 
  776. double tol;  /* the tolerance fit  */
  777. {
  778.     int i,
  779.         sym_op,                   /* A symmetry operation of the square */
  780.         ialpha,                   /* Intgerized scalling factor         */
  781.         ibeta,                    /* Intgerized offset                  */
  782.         best_ialpha,              /* best ialpha found                  */
  783.         best_ibeta,
  784.         best_sym_op,
  785.         best_domain_x,
  786.         best_domain_y, 
  787.         first_class, 
  788.         the_first_class, 
  789.         first_class_start,        /* loop beginning and ending values   */
  790.         first_class_end, 
  791.         second_class[2],
  792.         the_second_class, 
  793.         second_class_start,       /* loop beginning and ending values   */
  794.         second_class_end, 
  795.         range_sym_op[2],          /* the operations to bring square to  */
  796.         domain_sym_op;            /* cannonical position.               */
  797.  
  798.     struct classified_domain *node;  /* var for domains we scan through */
  799.  
  800.     double rms, best_rms,          /* rms value and min rms found so far */
  801.            sum=0, sum2=0;          /* sum and sum^2 of range pixels      */
  802.  
  803.  
  804.     /* keep breaking it down until we are small enough */
  805.     if (depth < min_part) {
  806.          quadtree(atx,aty, xsize/2, ysize/2, tol,depth+1);
  807.          quadtree(atx+xsize/2,aty, xsize/2, ysize/2,tol,depth+1);
  808.          quadtree(atx,aty+ysize/2, xsize/2, ysize/2,tol,depth+1);
  809.          quadtree(atx+xsize/2,aty+ysize/2, xsize/2, ysize/2,tol,depth+1);
  810.          return;
  811.     }
  812.  
  813.     /* now search for the best domain-range match and write it out */
  814.     best_rms = 10000000000.0;      /* just a big number */
  815.  
  816.     classify(atx, aty, xsize,ysize, 
  817.               &the_first_class, &the_second_class, 
  818.               &range_sym_op[0], &sum, &sum2, 1);
  819.  
  820.  
  821.     /* sort out how much to search based on -f and -F input flags */
  822.     if (fullclass_search) {
  823.          first_class_start = 0; 
  824.          first_class_end = 3;
  825.     } else {
  826.          first_class_start = the_first_class; 
  827.          first_class_end = the_first_class+1;
  828.     } 
  829.  
  830.     if (subclass_search) {
  831.          second_class_start = 0; 
  832.          second_class_end = 24;
  833.     } else {
  834.          second_class_start = the_second_class; 
  835.          second_class_end = the_second_class+1;
  836.     } 
  837.  
  838.     /* these for loops vary according to the optimization flags we set */
  839.     /* for subclass_search and fullclass_search==1, we search all the  */
  840.     /* domains (except not all rotations).                             */
  841.     for (first_class = first_class_start; 
  842.          first_class < first_class_end; ++first_class)
  843.     for (second_class[0] = second_class_start; 
  844.          second_class[0] < second_class_end; ++second_class[0]) {
  845.  
  846.        /* We must check each domain twice. Once for positive scaling,  */
  847.        /* once for negative scaling. Each has its own class and sym_op */
  848.        if (!only_positive) {
  849.           second_class[1] = 
  850.              class_transform[(first_class == 2 ? 1 : 0)][second_class[0]];
  851.           range_sym_op[1] = 
  852.              rot_transform[(the_first_class == 2 ? 1 : 0)][range_sym_op[0]];
  853.        } 
  854.  
  855.        /* only_positive is 0 or 1, so we may or not scan                */
  856.        for (i=0; i<(2-only_positive); ++i) 
  857.        for (node = the_domain[first_class][second_class[i]][depth-min_part]; 
  858.            node != NULL; 
  859.            node = node->next) {
  860.            domain_sym_op = node->the->sym;
  861.            /* The following if statement figures out how to transform  */
  862.            /* the domain onto the range, given that we know how to get */
  863.            /* each into cannonical position.                           */
  864.            if (((domain_sym_op>3 ? 4: 0) + (range_sym_op[i]>3 ? 4: 0))%8 == 0)
  865.              sym_op = (4 + domain_sym_op%4 - range_sym_op[i]%4)%4;
  866.            else 
  867.              sym_op = (4 + (domain_sym_op%4 + 3*(4-range_sym_op[i]%4))%4)%8;
  868.             
  869.            rms = compare(atx,aty, xsize, ysize,  depth, sum,sum2, 
  870.                                   node->the->dom_x, 
  871.                                   node->the->dom_y, 
  872.                                   sym_op, &ialpha,&ibeta); 
  873.     
  874.            if (rms < best_rms) {
  875.                    best_ialpha = ialpha;
  876.                    best_ibeta = ibeta;
  877.                    best_rms = rms;
  878.                    best_sym_op = sym_op;
  879.                    best_domain_x = node->the->dom_x;
  880.                    best_domain_y = node->the->dom_y;
  881.            }
  882.        }
  883.     }
  884.         
  885.     if (best_rms > tol && depth < max_part) {
  886.        /* We didn't find a good enough fit so quadtree down */
  887.        pack(1,(long)1,output);  /* This bit means we quadtree'd down */
  888.        quadtree(atx,aty, xsize/2, ysize/2, tol,depth+1);
  889.        quadtree(atx+xsize/2,aty, xsize/2, ysize/2, tol,depth+1);
  890.        quadtree(atx,aty+ysize/2, xsize/2, ysize/2, tol,depth+1);
  891.        quadtree(atx+xsize/2,aty+ysize/2, xsize/2, ysize/2, tol,depth+1);
  892.     } else {
  893.        /* The fit was good enough or we just can't get smaller ranges */
  894.        /* So write out the data */
  895.        if (depth < max_part)          /* if we are not at the smallest range */
  896.                pack(1,(long)0,output);/* then we must tell the decoder we    */
  897.                                       /* stopped quadtreeing                 */
  898.        pack(s_bits, (long)best_ialpha, output);
  899.        pack(o_bits, (long)best_ibeta, output);
  900.        /* When the scaling is zero, there is no need to store the domain */
  901.        if (best_ialpha != zero_ialpha) {
  902.           pack(3, (long)best_sym_op, output);
  903.           pack(bits_needed[depth-min_part], (long)(best_domain_y*
  904.             domain.no_h_domains[depth-min_part]+best_domain_x), output);
  905.       }
  906.      }
  907. }
  908.  
  909. /* ************************************************************* */
  910. /* Recursively partition an image, finding the largest contained */
  911. /* square and call the quadtree routine the encode that square.  */
  912. /* This enables us to encode rectangular image easily.           */
  913. /* ************************************************************* */
  914. partition_image(atx, aty, hsize,vsize, tol)
  915. int atx, aty, hsize,vsize;
  916. double tol;
  917. {
  918.    int x_exponent,    /* the largest power of 2 image size that fits */
  919.        y_exponent,    /* horizontally or vertically the rectangle.   */
  920.        exponent,      /* The actual size of image that's encoded.    */
  921.        size, 
  922.        depth; 
  923.    
  924.    x_exponent = (int)floor(log((double)hsize)/log(2.0));
  925.    y_exponent = (int)floor(log((double)vsize)/log(2.0));
  926.    
  927.    /* exponent is min of x_ and y_ exponent */
  928.    exponent = (x_exponent > y_exponent ? y_exponent : x_exponent);
  929.    size = 1<<exponent; 
  930.    depth = max_exponent - exponent;
  931.    quadtree(atx,aty,size,size,tol,depth);
  932.    if (size != hsize) 
  933.       partition_image(atx+size, aty, hsize-size,vsize, tol);
  934.  
  935.    if (size != vsize) 
  936.       partition_image(atx, aty+size, size,vsize-size, tol);
  937. }        
  938.   
  939. /* fatal is for when there is a fatal error... print a message and exit */
  940. void fatal(s)
  941. char *s;
  942. {
  943.      printf("%s\n",s);
  944.      exit(-1);
  945. }
  946.